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A phenomenological equation called Landau-Lifshitz-Baryakhtar (LLBar) equation, which could 
be viewed as the combination of Landau-Lifshitz (LL) equation and an extra “exchange damping” 
term, was derived by Baryakhtar using Onsager’s relations. We interpret the origin of this “exchange 
damping” as nonlocal damping by linking it to the spin current pumping. The LLBar equation is 
investigated numerically and analytically for the spin wave decay and domain wall motion. Our 
results show that the lifetime and propagation length of short-wavelength magnons in the presence 
of nonlocal damping could be much smaller than those given by LL equation. Furthermore, we find 
that both the domain wall mobility and the Walker breakdown field are strongly influenced by the 
nonlocal damping. 

PACS numbers: 75.78.Cd, 76.50.+g, 75.60.Ch 


I. INTRODUCTION 

The genuine complexity of magnetic and spintronic 
phenomena occurring in magnetic samples and devices 
imposes both fundamental and technical limits on the 
applicability of quantum-mechanical and atomistic the¬ 
ories to their modeling. To a certain degree, this chal¬ 
lenge can be circumvented by exploiting phenomenologi¬ 
cal theories based on the continuous medium approxima¬ 
tion. The theories operate with the magnetization (i.e. 
the magnetic moment density) and the effective magnetic 
field as generalized coordinates and forces respectively 
mm- The effective magnetic field is defined in terms 
of various magnetic material parameters, which are de¬ 
termined by fitting theoretical results to experimental 
data, and at least in principle, can be calculated using 
the quantum-mechanical or atomistic methods. However, 
solving the phenomenological models analytically is still 
a formidable task in the majority of practically impor¬ 
tant cases. The difficulty is primarily due to the pres¬ 
ence of the long range magneto-dipole interaction and 
associated non-uniformity of the ground state configura¬ 
tions of both the magnetization and effective magnetic 
field. Hence, the phenomenological models are solved in¬ 
stead numerically, using either finite-difference or finite- 
element methods realized in a number of micromagnetic 
solvers m- 

Traditionally, the software for such numerical micro- 
magnetic simulations of magnetization dynamics is based 
on solving the Landau-Lifshitz equation [IJ with a trans¬ 
verse magnetic relaxation term, either in the original 
(Landau) p] or ’’Gilbert” [8] form. Over time, dictated 
by the experimental and technological needs, the solvers 


have been modified to include finite temperature effects 
j9| and additional contributions to the magnetic energy 
(and therefore to effective magnetic field) ]TD) . The re¬ 
cent advances in spintronics and magnonics have led to 
the implementation of various spin transfer torque terms 
mm and periodic boundary conditions [T3fll5| . Fur¬ 
thermore, the progress in experimental investigations of 
ultrafast magnetization dynamics [lfij has exposed the 
need to account for the variation of the length of the mag¬ 
netization vector in response to excitation by femtosec¬ 
ond optical pulses, leading to inclusion of the longitudi¬ 
nal relaxation of the magnetization within the formalism 
of numerical micromagnetics m- Provided that a good 
agreement between the simulated and measured results 
is achieved, a microscopic (i.e. quantum-mechanical or 
atomistic) interpretation of the experiments can then be 
developed. 

The described strategy relies on the functional com¬ 
pleteness of the phenomenological model. For instance, a 
forceful use of incomplete equations to describe phenom¬ 
ena originating from terms missing from the model may 
result in false predictions and erroneous values of fitted 
parameters, and eventually in incorrect conclusions. The 
nature of the magnetic relaxation term and associated 
damping constants in the Landau-Lifshitz equation is of 
paramount importance both fundamentally and techni¬ 
cally. It is this term that is responsible for establish¬ 
ment of equilibrium both within the magnetic sub-system 
and with its environment (e.g. electron and phonon sub¬ 
systems), following perturbation by magnetic fields, spin 
currents, and/or optical pulses [IB] . Moreover, it is the 
same term that will eventually determine the energy effi¬ 
ciency of any emerging nano-magnetic devices, including 
both those for data storage [IS] and manipulation m- 
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In this report, we demonstrate how the phenomeno¬ 
logical magnetic relaxation term derived by Baryakhtar 
to explain the discrepancy between magnetic damping 
constants obtained from ferromagnetic resonance (FMR) 
and magnetic domain wall velocity measurements in di¬ 
electrics [20H22j can be applied to magnetic metallic sam¬ 
ples. We show that the Landau-Lifshitz equation with 
Baryakhtar relaxation term (Landau-Lifshitz-Baryakhtar 
or simply LLBar equation) contains the Landau-Lifshitz- 
Gilbert (LLG) equation as a special case, while also natu¬ 
rally including the contribution from the nonlocal damp¬ 
ing in the tensor form of Zhang and Zhang [23] and De 
Angeli [24| . The effects of the longitudinal relaxation 
and the anisotropic transverse relaxation on the magne¬ 
tization dynamics excited by optical and magnetic field 
pulses, respectively, in continuous films and magnetic ele¬ 
ments were discussed e.g. in Ref. Qj 5ESIES] ■ So, here we 
focus primarily on the manifestations of the Baryakhtar 
relaxation in problems specific for magnonics [T2 I and 
domain wall dynamics [27] [25]. This is achieved by in¬ 
corporating the LLBar equation within the code of the 
Object Oriented Micromagnetic Framework (OOMMF) 
[3], probably the most popular micromagnetic solver cur¬ 
rently available, and by comparing the results of simula¬ 
tions with those from simple analytical models. Specif¬ 
ically, we demonstrate that the Baryakhtar relaxation 
leads to increased damping of short wavelength spin 
waves and to modification of the domain wall mobility, 
the latter being also affected by the longitudinal relax¬ 
ation strength. 

The paper is organized as follows. In Sec. II, we re¬ 
view and interpret the Baryakhtar relaxation term. In 
Sec. Ill, we calculate and analyze the spin wave decay in 
a thin magnetic nanowire. In Sec. IV, we simulate the 
the suppression of standing spin waves in thin film. In 
Sec. V, we analyze the domain wall motion driven by the 
external field and compare the relative strength of con¬ 
tributions from the longitudinal and nonlocal damping. 
We conclude the discussion in Sec. VI. 


II. BASIC EQUATIONS 

In the most general case, the LLBar equation can be 
written as [201123] 

<9M 

— = —qM x H off + R (1) 

where q(> 0) is the gyromagnetic ratio and the relaxation 
term R is 

R = A r -H cff -A e , sp | 2]E ^. (2) 

l/J/ g *Xj Op 

Here and in the rest of the paper, the summation is 
automatically assumed for repeated indices. The two 
relaxation tensors A r and A e describe relativistic and 
exchange contributions, respectively, as originally intro¬ 
duced in Ref. [U. 


To facilitate comparison with the LLB equation as 
written in Ref. j 29j . the magnetic interaction energy of 
the sample is defined as 


w = w ll 


IM) (M 2 - M'e ) 2 

8x Ml 


( 3 ) 


where M e is the equilibrium magnitude of the magne¬ 
tization vector at a given temperature and zero micro- 
magnetic effective field, i.e. the effective field derived 
from the micromagnetic energy density w^, as used in 
standard simulations at constant temperature under con¬ 
dition |M| = M e = const (i.e. with only the trans¬ 
verse relaxation included). The second term in right- 
hand side of Eq. ^ describes the energy density induced 
by the small deviations of the magnetization length from 
its equilibrium value M e at the given temperature, i.e., 
|M 2 — Ml | M 2 , and \ is the longitudinal magnetic sus¬ 
ceptibility. Therefore, the associated effective magnetic 
field is 

h -' = -^ = h " + 4 (1 -" 2,m (4) 


where n = M/M e , H M is the effective magnetic field 
associated to uq,. Hereafter we assume that our system 
is in contact with the heat bath, so that the equilibrium 
temperature and associated value of M e and x remain 
constant irrespective of the magnetization dynamics. 

In accordance with the standard practice of both mi¬ 
cromagnetic simulations and analytical calculations, to 
solve LLBar equations ( l|4 ), one first needs to the cor¬ 
responding static equations obtained by setting the time 
derivatives to zero and thereby to derive the spatial dis¬ 
tribution of the magnetization in terms of both its length 
and direction. We note that, in general (e.g. as in the 
case of a domain wall), the resulting distribution of the 
longitudinal effective field and therefore also of the equi¬ 
librium magnetization length is nonuniform, so that the 
length is not generally equal to M e . With the static 
solution at hands, the dynamical problem is solved so 
as to find the temporal evolution of the magnetization 
length and direction following some sort of a perturba¬ 
tion. Crudely speaking, the effect of the relaxation terms 
is that, at each moment of time, the magnetization di¬ 
rection relaxes towards the instantaneous direction of the 
effective magnetic field, while the magnetization length 
relaxes towards the value prescribed by the instantaneous 
longitudinal effective magnetic field. The effective field 
itself varies with time, which makes the problem rather 
complex. However, this is the same kind of complexity as 
the one that has always been inherent to micromagnet¬ 
ics. The account of the longitudinal susceptibility within 
the LLBar equation only brings another degree of free¬ 
dom (the length of the magnetization) into the discus¬ 
sion. One should note however that the longitudinal sus¬ 
ceptibility has a rather small value at low temperature 
and so its account is only required at temperatures of the 
order of the Curie temperature. 
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We neglect throughout the paper any effects due to 
the anisotropy of relaxation, which could be associated 
e.g. with the crystalline structure of the magnetic ma¬ 
terial I2H1 E5|. This approximation is justified for poly¬ 
crystalline and amorphous soft magnetic metals, as has 
been confirmed by simulations presented in Ref. [25]. 
Hence, we represent the relaxation tensors as A r = A r I 
and A e = A e I where parameters A r and A e are the rela¬ 
tivistic and exchange relaxation damping constants and 
I is the unit tensor. Then, Eq. ([l]) is reduced to 

3 t M = - 7 M x H cff + A r H eff - A e V 2 H cff . (5) 

We separate the equations describing the dynamics and 
relaxation of the length and direction of the magnetiza¬ 
tion vector. Representing the latter as a product of its 
magnitude and directional unit vector M = Mm, we can 
write 


<9m dM 

M°— + m WT = - 7 M x H off + R. ( 6 ) 
at at 


We multiply this equation by m to obtain, 


dM 

~9t 


= m • R. 


( 7 ) 


Then, subtracting the product of equation |7| and m 
from equation Q, we obtain 


3m 

Ik 


— 7 m x H cff + — R_l 


( 8 ) 


where Rj^ = — m x (m x R). In the rest of the paper, we 
will use Ax = (A)x = A — (A ■ m)m to represent the 
component of the vector A that is perpendicular (trans¬ 
verse) to vector m. Note that only the perpendicular 
component of the torque contributes to cAm = dm/dt. 
For given temperature, M e is constant and we can define 
a = A r /( 7 M e ). In the limiting case of x —0, M —> M e 
and thus a is recognized as the Gilbert damping con¬ 
stant from the LLG equation. Let us now consider the 
case of A e ^ 0. The corresponding contribution to the 
relaxation term, which we denote here as Bgari can be 
written as 


B Bar = —A e V 2 H cff = -3J*, (9) 

where di = d/dxt and the quantity j, = — A e 3iH e ff has 
the form of some magnetization current density (magne¬ 
tization flux). 

For the following, it is useful to split the effective field 
into its perpendicular (relative to m) part (H^, “perpen¬ 
dicular field”) and parallel part (Hj| ff , “parallel field”), 

i.e., H off = and then to consider the associ¬ 

ated magnetic fluxes and torques separately. The mag¬ 
netic flux of j||,* = — A e 3jHj| ff and then the contribution 
of the associated torque T|| = — 3J|i j onto m is 

(t,|)x = -2A e 3,H|! ff 3 l m- A e H|J ff (V 2 m)x. (10) 


The perpendicular field can be represented as 


H 


_L _ 
eff 


l 

7 m 2 


M x 


3M 

Ik 


+ 0{ R) 


1 

7 


[m x 3 t m]. (11) 


So, we can write for the magnetization flux associated 
with the perpendicular field 


j_L,i = —(A e /7)3j(m x 3 t m). 


( 12 ) 


The right-hand side of Eq. (12) could be regarded as the 
torque generated by spin current pumping since mxftm 
can be considered as the exchange spin current [30], and 
then for the associated perpendicular torque t_l, we ob¬ 
tain, 


r± = -3jjx,i = -<jM e did t (m x 3 f m), (13) 

where we have introduced variable a = A e /( 7 M e ). We 
show that the torque (tj_)j_ could be written as (see Ap¬ 
pendix [A] for details) 

(tx)x = M e [m x (V • 3tm) — am x V 2 3tm] (14) 

where V is a 3 x 3 tensor E3ED, 

T> a p = 2cr(m x 3,;m) a (m x cAm)^ — a(dim ■ 3im)<5 a/ 3. 

(15) 

In the limit of x ~> 0, we assume Hjj ff = 0 and therefore 
obtain 


3 t m = — 7 m x H c ff — 7 am x (m x H c ff) 

-l-m x (V ■ 3*m) — urn x V 2 3f.m. (16) 

At the same time, Eq. ([ 8 ]) can then be written as 
c)yw 

— = - 7 m x Hoff - 7 m x (m x H® ff ), (17) 

where 

Hog- = aH eff - aV 2 H^ ff , (18) 


and H^j is the transverse component of the effective 
field. The first term in Eq. (18) is kept as Hoff since 
m x H c ff = m x H^. In practice, we use Eq. (17) rather 
than Eq. (16) for numerical implementation. As shown 
in Eq. (16) the damping terms contain both the form 
— m x V z O f m [Ml S21 and tensor form m x (V ■ 3 t m) [23]. 
Hence, we conclude that the exchange damping can 
be explained as the nonlocal damping, and Eq. @ is 
the phenomenological equation to describe the nonlocal 
damping. 

The intrinsic Gilbert damping is generally considered 
to have the relativistic origin T : ,'333- Phenomenologically, 
the Gilbert damping is local and the damping due to the 
nonuniform magnetization dynamics being ignored [ 8 ]. 
The exchange relaxation term in the LLBar equation de¬ 
scribes the nonlocal damping due to the nonuniform ef¬ 
fective field. Despite the complexity of various damping 
mechanisms, the spin current j in conducting ferromag- 
nets can be calculated, e.g. using the time-dependent 
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Pauli equation within the s-d model. The spin current 
is then given by j, = (g^_B?iGo/4e 2 )((9 t m x c^m), where 
Go is the conductivity [23], and thus the nonlocal damp¬ 
ing of the tensor form can be obtained m ini- As 
we can see from Appendix [A] this spin current densi¬ 
ties j, and j“ have the same form, and therefore, we 
can establish that er ~ ggs^G^/ 4e 2 M e . The spin cur¬ 
rent component $ ( se e Appendix |a| gives the term 
m x V 2 <9fm ES3, and the value of a can be therefore 
interpreted as a ~ (7 / goM e )(h/2) 2 n e T sc /m* , where n e 
is the conduction electron density, m* the effective mass 
and r sc is the transverse spin scattering time [34]. 

It is of interest to compare Eq. ([ 5 ]) with Landau- 
Lifshitz-Bloch (LLB) equation 29], which could be writ¬ 
ten as 

dU tt , 7 “H r tt 1 „ s 

— = -qn x H c ff H-^ n • H cff n -—n x (n x H eff ) 

at n z n z 

(19) 

where n = M /M e (T) is the reduced magnetization and 
M e (T) is the equilibrium magnetization value at tem¬ 
perature T. The effective field H c ff contains the usual 
micromagnetic contributions Hj nt as well as the contri¬ 
bution from the temperature, 

H eff = H int + ^(l-n 2 )n (20) 

2X|| 


where m e = M e (T)/M e ( 0) and xy = dm/dH with m = 
M/M e ( 0) [29,. By substituting Eq. fl20| into Eq. (19), 
one arrives at 

<9n 

— = -qn x H int + qay (H int )|| + qau(H int )_L 

+ ^(l^ 2 )n. (21) 

2 X|| 

Meanwhile, if we neglect the A e term in Eq. |5]) and insert 
the effective field Eq. into Eq. ([ 5 ]), we obtain 

dw. A 

= -qn X H int + qArHint + (! - n 2 ) n. (22) 


As we can see, Eq. (22) is a special case of LLB equation 
with the assumption that a± = a\\ = A r /(qM e ) and \ = 
M e (0)%||. However, the LLB equation does not contain 
the A e -term (nonlocal damping term) which is the main 
focus in this work. 


III. SPIN WAVE DECAY 


size is 1 x 2 x 2nm 3 . The magnetization aligns along 
the e x direction for the equilibrium state and the pa¬ 
rameters are typical of Permalloy: the exchange con¬ 
stant A = 1.3 x 10 -11 J/m, the saturation magnetization 
M e = 8.6 x 10 5 A/m and the Gilbert damping damping 
coefficient a = 0.01. The spin waves are excited locally 
in the region 0 < x < 2 nm, and to prevent the spin wave 
reflection the damping coefficient is increased linearly [35] 
from 0.01 at x = 1802 nm to 0.5 at x = 2002 nm. 



x (nm) 


FIG. 1. The spin wave amplitude decay along the rod, for a 
spin wave was excited locally by applying a microwave H = 
Ho sin(27r ft)e y of frequency / = 30 GHz and amplitude Ho = 
1000 Oe in the region 0 < x < 2 nm. The data were htted 
using Eq. ( |23| ) with (3 = 0.02 and a = 0.01. 


Figure [I] illustrates the spin wave amplitude decay 
along the rod. The y component of magnetization unit 
vector uiy data for 30 < x < 1800 nm were fitted using 
(23) to extract the wave vector k and the decay constant 


A, and good agreement is observed due to the effective 
absence of spin wave reflection. We use data after having 
computed the time development of the magnetization for 
4 ns to reach a steady state. The injected spin wave en¬ 
ergy is absorbed efficiently enough within the right 200 
nm of the rod due to the increased damping. 

To analyze the simulation data, we exploit the uniform 
plane wave assumption with its exponential amplitude 
decay due to energy dissipation, i.e. magnetization with 
the form e -\x ^ w j lere \ j s tl le characteristic pa¬ 

rameter of the spin wave damping. For a small amplitude 
spin wave propagation we have [551 


m = e x 


_ mo e i(kx-ut) e -\x 


(23) 


where |mo| AC 1, and the effective field of the long rod 
can be expressed as 

H c ff = H s m x e x + HV 2 m, (24) 


To perform the micromagnetic simulation for the spin 
wave decay, we have implemented Eq. <0 as an ex¬ 
tension for the finite difference micromagnetics package 
OOMMF. A new variable /3 for the exchange damping is 
introduced with a = /3G, where G is a coefficient to scale 
0 to the same order as a. In practice, G was chosen to 
be G = A/^ioM 2 ). 

The simulation geometry has dimensions L x = 
2002 nm, L y = 2 nm and L z = 2 nm, and the cell 


where the ‘easy axis’ anisotropy field H s m x e x originates 
from the demagnetizing field, and the constant D mea¬ 
sures the strength of the exchange field, 


H a 


2 K 
go M e 



2 A 

y 0 M e ' 


(25) 


To test the spin wave decay for this system, a sinusoidal 
field H = H 0 sin(27r ft)e y was applied to the rod in the 
region 0 < x < 2 nm to generate spin waves. 
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is plotted using Eq. ([29]) with (3 = 0.01 and a = 0.01, 
which shows a good approximation for the simulation 
data. Besides, this exchange damping could be impor¬ 
tant in determining the nonadiabatic spin torque. We 
could establish the value of (3 using the existing exper¬ 
imental data, such as the transverse spin current data 
[31] gives j3 ~ 0.1 which hints the lifetime and propaga¬ 
tion length of short-wavelength magnons could be much 
shorter than those given by the LLG equation EU- 


IV. SUPPRESSION OF STANDING SPIN 
WAVES 


FIG. 2. The spin wave decay constant-wave vector product 
Afc as a function of the frequency for different (3 values. The 
slateblue line was drawn using Eq. (|29|) for the case /3 = 0.01. 


Figure [2] shows the product of spin wave decay con¬ 
stant A and wave vector A; as a function of the frequency. 
The dependence is linear for the f3 = 0 case, which is 
agreement with the zero adiabatic spin torque case |36| . 
The addition of a nonzero (3 term leads to a nonlinear 
relation, and the amplitude of the spin wave decay con¬ 
stant that is significantly larger than that given by the 
linear dependence. We also performed the simulation for 
X > 0 case by using Eq. <®> which shows that (3 term is 
the leading factor for this nonlinearity (the relative error 
is less than 1% for x = 1 x 10 -3 ). To analyze the nonlin¬ 
ear dependence, we introduce the complex wave vector 
k, 


k = k + A i. (26) 

By linearizing Eq. <H3> and setting the determinant of the 
matrix to zero we obtained (see Appendix |B| for details): 

(w + cDq -f- iCj i)(w — Wo + u hi) = 0, (27) 


where wo = "f(H s + Dk 2 ) and uq = awo + (3Gk 2 ujQ. The 
second term of the Eq. ( [27| is expected to be equal to 
zero, i.e., oq — ico + iuj o = 0. There are two scenarios to 
consider: First is the (3 = 0 case, kX could be extracted 
by taking the imaginary part of k 2 at Eq. (26): 


kX = - Im 


{*} 


au 


2(1 + a 2 )~/D' 


(28) 


The linear dependence of kX as a function of frequency 
matches the data plotted in Fig. [2] For the (3 > 0 case, 
solving Eq. (27) yields in the linear with respect to the 


damping constants approximation, 


kX 


OJ 

2r/D 


(a + (3Gk 2 ) 


(29) 


where the dispersion relation for the rod is w = 7 (H s + 
Dk 2 ). Eq. (29) shows there is an extra k 2 term associated 
with the exchange damping term besides the linear de¬ 
pendence between kX and w. The slateblue line in Fig. [2] 


In the presence of nonlocal damping, the high fre¬ 
quency standing spin waves in the thin films are sup¬ 
pressed m If the magnetization at the surfaces are 
pinned, the spin wave resonance can be excited by a uni¬ 
form alternating magnetic field [38]. With given out-of- 
plane external held H z in the ^-direction, the frequencies 
of the excited spin waves of the him are given by [55], 

\ 2 

( 30 ) 

where d is the him thickness, u>o = ”/(H z — M e ), = 
yMe and A ex = y/2A/(fioM£). The excited spin wave 
modes are labeled by the integer n, and the odd n has 
a nonvanishing interaction with the given uniform alter¬ 
nating magnetic held [55], 

300 

» 200 
a 100 

0 


2 4 6 8 10 12 

Frequency f (GHz) 



FIG. 3. Imaginary parts of the dynamical susceptibility \vv 
of the him for (a) thickness d = 300 nm, and (b) thickness 
d = 60 nm. 



To reduce the simulation time, we consider a system 
with cross-sectional area 4x4 nm 2 in xy-plane and apply 
the two-dimensional periodic boundary conditions [14] to 
the system. We use the Permalloy as the simulation ma¬ 
terial with external held H z = 1 x 10 6 A/m and the cell 
size is 4 x 4 x 2 nm 3 . Instead of applying microwaves to 
the system, we calculate the magnetic absorption spec¬ 
trum of the him by applying a sinc-function held pulse 
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h = /losinc(woi) to the system [50]. With the collected 
average magnetization data, the dynamic susceptibility 
X is computed using Fourier transformation m- For ex¬ 
ample, the component Xyy is computed using rn y when 
the pulse is parallel to the y- axis. 

Figure [3^ a) shows the imaginary part of the dynamic 
susceptibility Xyy f° r a film with d = 300 nrn. As we 
can see, the spin wave of modes n = 1,3, 5,... are ex¬ 
cited, and the influence of the “exchange damping” is 
small. However, the presence of the “exchange damp¬ 
ing” suppresses the spin wave excitation (n > 1 mode) 
significantly for the film with thickness d = 60 ran, as 
shown in Fig. [3])b). The reason is because the damping 
of the standing spin waves is proportional to fc 4 in the 
presence of exchange damping EZj. 


V. DOMAIN WALL MOTION 


We implemented Eq. ([5]) in a finite element based mi- 
cromagnetic framework to study the effect of parallel re¬ 
laxation process on domain wall motion. The simulated 
system for the domain wall motion is a one-dimensional 
(ID) mesh with length 20000 nm and discretization size 
4nm, a head-to-head domain wall is initialized with its 
center near x = 500 nm. In this section, the demagnetiz¬ 
ing fields are simplified as Hu = — AM and the demag¬ 
netizing factors are chosen to be N x = 0, N y = 0.2 and 
N z = 0.8, respectively. The domain wall moves under 
the applied field for 50 ns and the domain wall velocities 
at different external field strengths are computed. Fig¬ 
ure H] shows the simulation results of domain wall motion 
under external fields for different susceptibilities with¬ 
out consideration of exchange damping, i.e., /3 = 0. For 
Nickel and Permalloy, the longitudinal susceptibility is 
around 10~ 4 at room temperature and increases with 
the temperature up until the Curie point [25] • We find 
that the longitudinal susceptibilities have no influence on 
the maximum velocity but change the Walker breakdown 
field H w significantly. The domain wall velocity in the 
limit x —^ 0 is almost the same with the case of x = 10 —4 , 
which could be explained by the relation that t he d iffer- 
ence proportional to the ratio of ( x/ a ) 2 i n Eq. (53). 


To investigate the effect of longitudinal magnetic sus¬ 
ceptibility x an d exchange relaxation damping a on the 
domain wall motion, we use the remainder of this sec¬ 
tion for analytical studies. We start from the constant 
saturation magnetization of one-dimensional domain wall 
model, such as the ID head-to-head wall [52]. The static 
ID domain wall profile can be expressed as 


= — tanh 


x — q 

A 


rrit = sech 


x — q 

A 


(31) 


where rrit is the perpendicular component of the unit 
magnetization vector, A is the wall width parameter and 
q is the position of the domain wall center. 

We consider the case that the system is characterized 
by two anisotropies, easy uniaxial anisotropy K and hard 



FIG. 4. Simulations results of domain wall velocities for var¬ 
ious susceptibilities. The parameters used are: a = 0.001, 
P = 0, N v = 0.2 and N z = 0.8. The vertical dash lines are 
the breakdown fields computed using Eq. (531. 


plane anisotropy K± , which originate from demagnetiza¬ 
tion. The aim is to analyze the impact of the longitudinal 
magnetic susceptibility under the ID domain wall model, 
the demagnetization energy density could be written as 


E» 


= M 

m? x m r 


K± 71 r‘2 


(32) 


where I\ = (l/2)(N y — N x )y 0 Ml and K± = (1/2 )(N Z — 
N y )y 0 M 2 . In the limit case x -t 0 case, the effective 
anisotropy energy density E an can be rewritten as 

E' an = Ksin 2 9(1 + Ksin 2 ip), (33) 


where m = (cos 9, sin 6 cosy), sin 6 sin <p) is used and 
n = K±/K is the ratio between hard plane anisotropy 
strength and easy uniaxial anisotropy strength. 

The dynamics of the domain wall with ID profile can 
be described using 3 parameters |43j : the domain width 
A, the domain wall position q and the domain wall tilt 
angle </>. In this domain wall model, one can assume 
that tp(x, t) = <p(t) is only a function of time. Thus, the 
magnetization profile for the head-to-head domain wall 
is given by 


9(x,t) = 2ta,n 1 exp(^^y^j , <p(x, t) = 

( 34 ) 

Using the magnetization unit vector to calculate the 
exchange energy is a good approximation for the case 
X<1, thus, the total energy density can be rewritten as 


Etot 


y 0 ( M 2 - M 2 e f 
8 X Ml 


m), 


where 


/ \ , vO o i^_L o 

] = M2 (Vm) ~ TW m * + 772 m *- 


Ml 


Ml 


(35) 


( 36 ) 
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Within the ID domain wall profile, H m , the longitudinal 
component of the effective field is obtained: 

H m = m • H eff = ^2 ( M 2 e - Ad 2 ) - 2AdP sin 2 9 (37) 


where P is defined as 

1 


P = 


Mo Ml [A 2 


A 


+ K(1 + k sin 2 <j>) 


(38) 


As we can see, P is a function of the tilt angle <j> and the 
domain wall width A. At the static state, H m should 
equal zero, i.e., dM/dt = 0, which gives 

Ad 2 = (1 - 4yPsin 2 9)Ad 2 . (39) 


Eq. (39) shows that the difference between magnetization 
length Ad and Ad e reaches its maximum at the center of 
the domain wall due to the effect of the exchange field, 
which also peaks in the centre of the domain wall. Ac¬ 
cording to Eq. p9| , we can estimate that the magneti¬ 
zation length difference is SAd m — 2%Psin 2 9 for \ 1 

case. Figure [5] shows the magnetization length differ¬ 
ences of a ID domain wall for various \i if can be seen 
that this approximation for 5Ad agrees very well with the 
simulation results. 



FIG. 5. Simulation results of the magnetization length differ¬ 
ence 5M for a ID domain wall located at x = 500 nm with 
M e = 8.6 x 10 5 A/m and A = 1.3 x 1CP 11 J/m. The demagne¬ 
tizing factors are selected to be N x = 0 and N v = N z =0.5. 


In the dynamic case, H m is not equal to zero. If we 
wrote Eq. (371 as H m = FAd, we can find that the non¬ 
trivial term that contributes to H m is 


F = -^(1 - Ad 2 /Ad 2 ) - 2P sin 2 9. (40) 

As an approximation for H m , we expect dF/dt = 0 [44j . 
which gives 


tt AP x q 2 

H m = -T- m t m x 

A a 7 


(41) 


by the domain wall velocity q only. We employ the La- 
grangian equation combined with dissipation function F 
to compute the domain wall dynamics m- The Lagrange 
equations are 


ac _ _d ddC\ 
dx~ dt \djc) 

where X refers to q, </> and A. 
is defined by F = f F dx where 



(42) 


The dissipation function 


F = i Mo M e7 [aH 2 ff + ff(VH c j) 2 ] (43) 


is the dissipation density function. 


A. Parallel relaxation 


We neglect the exchange damping term with assump¬ 
tion that a <C aA 2 and arrive at 

F = l -a^Ad el H 2 eS = iaMoM e7 (Hi + H 2 m ). (44) 

where Hj^ and H,„ are the perpendicular and parallel 
components of the effective field. If we also assume that 
a ~ x < !- can be approximated by Eq. (|ll|), 


Hi 


= -4 m 2 = \( 9 2 

' 1 * rsjA 


sin 


(45) 


Substituting Eq. (411 and Eq. (|45[) into Eq. M4| and 


integrating over space, we obtain 


F = 


a no Ad e 


i 2 A + J(l + Q) 


(46) 


A = 


where we have ignored the A term. This term leads to 
the optimal domain wall width m- 

\Ja/(K + AAsin 2 (/)) (47) 

and for k = 0 the optimal domain wall width reduces 
to Ao = \J’A/K. In what follows, the domain wall 
width parameter A(t) is approximated by the optimal 
wall width. The parameter P is then given by 


P = 


2 /v (1 + « sin 2 (j>) 


2 A 
HoAdl A2’ 


(48) 


HoAd 2 

and it is straightforward to find its minimum Pq = 
2K/(hqM 2 ), which corresponds to A = A 0 . 

The introduced paramter Q in Eq. (46) is given by 
Q = (32/15 )P 2 (x/ a ) 2 an d its value is determined by the 
ratio of x an d a, which could be ~ 1 although we assume 
X ~ a 1. Following the treatment of Ref. m , the 
integrated Lagrangian action C is given by 

Ho M e 


C = J (E t , 
2 A 


7 


-(/cos9) dx 


= ^ + 2A/v (1 + ft sin 2 (/>)(! - V) 


( 49 ) 


^ 2 HoM e dd a q 


2 HoM e 


Hi 


In this approximation, we have ignored the terms con¬ 
taining dP/dt and thus the amplitude of dd m is influenced 
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where iM}M e (j> cos 0/7 is the Berry phase term [15], V = 
8yP/3 is a result of the varying magnetization that in¬ 
troduced a pinning potential. However, the potential is 
fairly small and therefore is negligible since V <C Q. By 
substituting Eq. (|49|) and Eq. (44) into Eq. (|42|) , 


</> + a— (1 + Q) — 7 H a , 

Q ; H k . 

—— a® = 7 — sin 2<p. 
A Y ' 2 Y 


(50) 


where H k = 2K±/(noM e ). The domain wall dynamics 
is governed by Eq. (501, by eliminating q we obtain an 
equation about <f >, 


7 


1 + a 2 ( 1 + <2) 


[Pa — H w (l + Q) sin 2</>] (51) 


where H w = aH k /2 is the Walker breakdown field. From 
Eq. (51) we can find that the critical value of is ap¬ 


proximately equal to 7 r /4 if Q <C 1, which leads to the 
maximum value of P to be P\ = 2K(1 + k/2)/(^j, 0 M 2 ). 
There exists an equilibrium state <j >* such that = 0 if 
H a <H W (1+Q), 


sin( 2 </>*) = h = 


Bn 


H w (l + Q) ’ 


(52) 


which means the Walker breakdown field H' w for % > 0 
case is increased to H w ( 1 + max{Q}), i.e., 


HL = H„ 


1 + —Pf - 


15 




(53) 


where P\ is the maximum vlaue of P. For this steady- 
state wall motion, the domain wall velocity is 


lH a A* 
a 1 + Q(A*) ’ 


(54) 


where 


A* 


A 0 /^l + ^(l-vT^). 


(55) 


Therefore, A* -7 Aq in the limit case H a — > 0, and the 
domain wall mobility /r is given by 


M = 


7 A 


a 


1 + - P 2 (-) 2 
1+ 15^°V 


(56) 


where Po is the minimum value of P. In Fig. [4] the cor¬ 
responding Walker breakdown fields are plotted in ver¬ 
tical dash lines, which gives a good approximation for 
% = 5 x ICE 4 and % = 1 x 10 ~ 4 cases. The simulation 
results show that the Walker breakdown field H w could 
be changed significantly if the longitudinal susceptibility 
is comparable to the damping constant. 



Applied field (A/m) 


FIG. 6. Simulations results of domain wall velocities for the 
limit case that y —¥ 0 with various exchange dampings. The 
parameters used are: a = 0.005, N y = 0.4 and N z = 0.6. The 
vertical dash lines are the breakdown fields computed with 
Eq. {60 1. 


B. Nonlocal damping 


In this part we consider the domain wall motion in¬ 
fluenced by exchange damping for the case that y —> 0 . 
The dissipation density function (43) thus becomes 

P = ^ 0 M e7 [aHi + a{VH e ) 2 + a(VP^) 2 ] (57) 

where Hg and H^ are the two components of the effec¬ 
tive field, and Hj^ is computed using Eq. (451. After 
calculation we obtain 


T = 


HoM e 

7 


; 9 , . 1 <7 , q 2 , 1 (J , 

^’ (aa+ 3d )+ i (Q+ 3Ai ) 


(58) 


We take the same Lagrangian action (49) for y = 0 and 
arrive at 


^+ ( a+^)f =7#a, 

q , . cr • H k . 

A' (a+ 3A = 


(59) 


Similarly, the corresponding Walker breakdown field 
changes to 


1 


H'w ~ 7 ;Hk ( a + 3 ^2 


1 a 


(60) 


where Ai = A 0 \/l/(l + re/2). The domain wall mobility 
is given by 


1 If la_\ 

M 7 A 0 V + 3 A o/ 


(61) 


As we can see, the nonlocal damping term a influences 
the domain wall motion as well, and we can establish 
that a/A' 2 = /?(1 + n/2)K/(n 0 M 2 ) oc j3. Therefore, for 
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the scenarios that K ~ fi 0 M 2 , the contributions from 
the Gilbert and nonlocal damping are of the order of 
magnitude for both the domain wall mobility and Walker 
breakdown field. 


Figure [ 6 ] shows the domain wall velocities for domain 
wall motion driven by external fields in the limiting case 
of x —>■ 0. The simulation results are based on a one¬ 
dimensional mesh with length 10000 nm with cell size of 2 
nm. The damping a is set to 0.005 and the demagnetizing 
factors are chosen to be N x = 0, N y = 0.4 and N z = 0.6. 


As predicted by Eq. (60), the nonlocal damping /3 leads to 
an increment of the Walker breakdown field, and Eq. (60) 
fits the simulation results very well. 


VI. SUMMARY 


We explain the “exchange damping” in the Landau- 
Lifshitz-Baryakhtar (LLBar) equation as nonlocal damp¬ 
ing by linking it to the spin current pumping, and there¬ 
fore the LLBar © can be considered as a phenomeno¬ 
logical equation to describe the nonlocal damping. In the 
presence of nonlocal damping, the lifetime and propaga¬ 
tion length of short-wavelength magnons could be much 
shorter than those given by the LLG equation. Our simu¬ 
lation results show that the spin wave amplitude decays 
much faster in the presence of nonlocal damping when 
spin waves propagate along a single rod. The analyti¬ 
cal result shows that there is extra nonlinear dependence 
scaling with k 2 between A k (the product of spin wave 
decay constant A and wave vector k) and frequency w 
due to the nonlocal damping. Using the micromagnetic 
simulation based on the LLBar equation, we show that 
the difference between magnetization length M and M e 
reaches its maximum at the center of the domain wall. 
For the cases that x ~ 01 where x is the longitudinal 
magnetic susceptibility and a is the Gilbert damping, the 
Walker breakdown field will increase significantly. By us¬ 
ing a ID domain wall model, we also show that both the 
domain wall mobility and the Walker breakdown field are 
strongly influenced by the nonlocal damping as well. 
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Appendix A: Derivation of equation (16) 


We split the perpendicular spin current jj_,j into two 
components, 

j_L,i = j“ +ji 5 (Al) 

where we write A e /7 as d, 

j“ = ~a(d t m x d t m) (A2) 

Ji = —d(m x did t m) (A3) 

The torque r a generated by spin current j“ is given by 
To. = {diji)±, Le -> 

T a = dm x [cj,rn x (d t m x c^m)] (A4) 

where we have used the identities m-c^rn = — dixn. d t va 
and m-didim = — Meanwhile, the correspond¬ 

ing torque 17 can be computed by = (9jji)j_, which 
gives 

Tf, = r a — a(di m ■ <9;m)m x — dm x V 2 9tm (A5) 

Note that r a = d<9,m[((9 t m x <9,m) • m] can be changed 
into the tensor form, 


Ta = m x (X>° • d t m), 


(A 6 ) 

where 



V° a p = d(m x <9 4 m) Q (m x c 

^m)^. 

(A7) 

Therefore, we obtain for r„ + 17 , 



T a + T b = m x (P ' 9 t m) — dm x 

<1 

to 

3 

(A 8 ) 


where V is a 3 x 3 tensor, 

V a p = 2 d(m x cAm) a (m x c^m)^ — d(dim ■ diTt\)8 a p. 

(A9) 


Appendix B: Derivation of equation (27) 


We introd uce a new variable s to represent the second 
term in the (23), i.e., s = m 0 e l( - fex_ “ Jt \ so we have 


m = e x + s, 


(Bl) 
(B2) 

H cff = H s ( 1 + 4)e x - DPs (B3) 

where s' x « (l/2)(s^, — s 2 ). Considering the fact |s| <C 1 
and neglect the high order term s 2 , one obtains = 
— (H s + Dk 2 )s and thus 

d s 


H b 


eff 


— C G x 


where 


c = aH s ( 1 + s' x ), 
d = —f3Gk 2 (Dk 2 + H s ) - aDk 2 


(B4) 

(B5) 

(B 6 ) 
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Substituting the above equations into we have where / = H s (l+s' x ) + Dk 2 . Neglecting high order terms 

such as s 2 and s x s y we obtained, 


7 [aH s — d) — ito u>o 

Sy 


O' 

—Wq 7 (aH s — d) — iw 

_s z 


0 


(B8) 


iuj 

Sx 

Sy 

= f 

1 

CO r-j 

+ (c — d) 

'-(S 2 y+S 2 z)' 
(1 + s x) s y 

Therefore, Eq. (271 can be obtained by setting the deter 
, (B7) minant of the matrix in (B8) to zero. 

7 

Sz 


Sy 


(1 + s x )s z 
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